function [p] = solve_contraction(fun,p0)

diff = 10;
tol = 10^-7;
p = p0;
t = 0;

check_contraction = 1;

while diff>tol && check_contraction==1

    p_new = fun(p);
    diff_new = max(abs(p-p_new));
    
    check_contraction = (diff_new < diff);
    
    p = p_new;
    diff = diff_new;
    t = t+1;

end

if check_contraction==0
    
    p = NaN(size(p0));
    
end

end

